Not really random but definitely confusing

---

I. “Random” numbers in R

In R random numbers are not truly random. Instead, they are generated using a seed value (set.seed) and a pseudo-random number generator (the default in R is Mersenne-Twister, but there are others). This RNG is a function that uses the seed to generate a deterministic sequence of numbers that approximates a sequence of truly random numbers, which is stored in a vector called .Random.seed. The first value of .Random.seed identifies the RNG. The length and structure of the remaining portion of the vector depends on the RNG itself. For Mersenne-Twister, the total vector is 625 elements long, the second of which identifies how many of the 623 values have been “used up” in creating output from functions that call .Random.seed (runif, rnorm, sample, etc). Once the end of the vector is reached, a fresh set of 624 values replace the first set and the position counter is set to 1. Oddly the initial position counter starts on the 623rd value of the RNG output when first created with set.seed. Invoke set.seed again with a new seed value and you will get an entirely new .Random.seed starting at position 624 (the end of the first set). Invoke set.seed with the same seed value and you will recreate .Random.seed from the beginning.

One tricky bit is how RNG’s work in regards to parallelized R workers. If you fork your master process to create R workers then all the workers will have the same RNG and seed as the R master process, which can be bad if you want the workers to be independent with regards to random number generation. To avoid this skip the fork option (see below) and each worker will default to the Mersenne-Twister RNG and a seed that depends on the system time and the worker’s process identifier (or PID). This insures that each worker’s .Random.seed vector is unique (but also non-reproducible since it is based on the time and PID). The upshot for running JAGS models is that you can now generate random initial values within each worker and know that they are going to be different from initial values generated in another worker.

set.seed(1)

# BAD!

cl <- parallel::makeCluster(3, type = "FORK")

out <- clusterEvalQ(cl, {
  a <- runif(3)
  x <- head(.Random.seed)
  return(x)

})

out

## [[1]]
## [1]         403           3  1654269195 -1877109783  -961256264  1403523942
## 
## [[2]]
## [1]         403           3  1654269195 -1877109783  -961256264  1403523942
## 
## [[3]]
## [1]         403           3  1654269195 -1877109783  -961256264  1403523942

stopCluster(cl)

# GOOD!

cl <- parallel::makeCluster(3)

out <- clusterEvalQ(cl, {
  a <- runif(3)
  x <- head(.Random.seed)
  return(x)

})

out

## [[1]]
## [1]         403           3 -1225896594    25615129  1820565419  -544883302
## 
## [[2]]
## [1]         403           3   586075502   535583513 -1083826261  -132563558
## 
## [[3]]
## [1]         403           3  -722416274 -1126540519   548937643   624377242

stopCluster(cl)

II. Simulate and fit some data

OK, so we quickly simulate some data for a random slopes random intercepts model.

N <- 10
SITES <- 20

meth <- meth.lmer <- c()

for (i in 1:SITES) {

  temp <- rbeta(N, 1, 1)
  meth <- cbind(meth, temp)
  meth.lmer <- c(meth.lmer, temp)

}

ages <- abs(round(rnorm(N, 30, 10)))

ages.lmer <- rep(ages, SITES)
sites.lmer <- sort(rep(seq(1:SITES), N))

AGES <- matrix(rep(ages, SITES), nrow = N, byrow = FALSE)

M_data <- list(
  y = AGES,
  x = meth,
  L = nrow(meth),
  S = ncol(meth))

We then recover parameter values with lmer. We will use these values to pick some initial values.

M1 <- lmer(ages.lmer ~ meth.lmer + (1 + meth.lmer | sites.lmer))

summary(M1)

## Linear mixed model fit by REML ['lmerMod']
## Formula: ages.lmer ~ meth.lmer + (1 + meth.lmer | sites.lmer)
## 
## REML criterion at convergence: 1567.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4485 -0.4394 -0.0129  0.7966  1.7574 
## 
## Random effects:
##  Groups     Name        Variance  Std.Dev.  Corr
##  sites.lmer (Intercept) 0.000e+00 0.000e+00     
##             meth.lmer   3.505e-15 5.921e-08  NaN
##  Residual               1.540e+02 1.241e+01     
## Number of obs: 200, groups:  sites.lmer, 20
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   37.655      1.798  20.947
## meth.lmer     -4.761      3.172  -1.501
## 
## Correlation of Fixed Effects:
##           (Intr)
## meth.lmer -0.873

coef(M1)

## $sites.lmer
##    (Intercept) meth.lmer
## 1      37.6548 -4.760572
## 2      37.6548 -4.760572
## 3      37.6548 -4.760572
## 4      37.6548 -4.760572
## 5      37.6548 -4.760572
## 6      37.6548 -4.760572
## 7      37.6548 -4.760572
## 8      37.6548 -4.760572
## 9      37.6548 -4.760572
## 10     37.6548 -4.760572
## 11     37.6548 -4.760572
## 12     37.6548 -4.760572
## 13     37.6548 -4.760572
## 14     37.6548 -4.760572
## 15     37.6548 -4.760572
## 16     37.6548 -4.760572
## 17     37.6548 -4.760572
## 18     37.6548 -4.760572
## 19     37.6548 -4.760572
## 20     37.6548 -4.760572
## 
## attr(,"class")
## [1] "coef.mer"

III. JAGS code

{
sink('Meth_1.jags')

cat("
    
model {
    
# y <- AGES (matrix - i x j)
# x <- METH (matix - i x j [each row is different individual, each column different site]
# hierarchical model - age ~ methylation at each site
# site-level alphas share a distribution, site-level betas share a distribution

#i is individual
#j is site
#L is number of individuals
#S is number of CPG sites

# priors

for (j in 1:S) {

  alpha[j] ~ dnorm(mu_a, tau[2])
  beta[j] ~ dnorm(mu_b, tau[3])

}
    
mu_a ~ dnorm(0, 0.001)
mu_b ~ dnorm(0, 0.001)

for (i in 1:3){

  sigma[i] ~ dunif(0, 1000)
  tau[i] <- pow(sigma[i], -2)    

}

# likelihood
    
for (i in 1:L) {
  for (j in 1:S) {
    
    y[i, j] ~ dnorm(mu[i, j], tau[1])
    mu[i, j] <- alpha[j] + beta[j] * x[i, j]

  }
}
    
}", fill = TRUE)
sink()
}

IV. Sequential chains with different initial values

We select different initial values for each chain fit the model in JAGS.

sigma <- NA

sigma[1] <- summary(M1)[[11]]
sigma[2] <- attributes(summary(M1)$varcor$sites.lmer)$stddev[1]
sigma[3] <- attributes(summary(M1)$varcor$sites.lmer)$stddev[2]

inits1 <- list(
  alpha = coef(M1)$sites.lmer[,1], 
  beta = coef(M1)$sites.lmer[,2],
  mu_a = summary(M1)[[10]][1, 1], 
  mu_b = summary(M1)[[10]][2, 1], 
  sigma = sigma + 1)
    
inits2 <- list(
  alpha = coef(M1)$sites.lmer[,1] + 5, 
  beta = coef(M1)$sites.lmer[,2] + 5,
  mu_a = summary(M1)[[10]][1, 1] + 5, 
  mu_b = summary(M1)[[10]][2, 1] + 5, 
  sigma = sigma + 2)

inits3 <- list(
  alpha = coef(M1)$sites.lmer[,1] + 10, 
  beta = coef(M1)$sites.lmer[,2] + 10,
  mu_a = summary(M1)[[10]][1, 1] + 10, 
  mu_b = summary(M1)[[10]][2, 1] + 10, 
  sigma = sigma + 3)

(F_Inits <- list(inits1, inits2, inits3))

## [[1]]
## [[1]]$alpha
##  [1] 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548
##  [9] 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548
## [17] 37.6548 37.6548 37.6548 37.6548
## 
## [[1]]$beta
##  [1] -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572
##  [8] -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572
## [15] -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572
## 
## [[1]]$mu_a
## [1] 37.6548
## 
## [[1]]$mu_b
## [1] -4.760572
## 
## [[1]]$sigma
## [1] 13.41029  1.00000  1.00000
## 
## 
## [[2]]
## [[2]]$alpha
##  [1] 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548
##  [9] 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548
## [17] 42.6548 42.6548 42.6548 42.6548
## 
## [[2]]$beta
##  [1] 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279
##  [8] 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279
## [15] 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279
## 
## [[2]]$mu_a
## [1] 42.6548
## 
## [[2]]$mu_b
## [1] 0.2394279
## 
## [[2]]$sigma
## [1] 14.41029  2.00000  2.00000
## 
## 
## [[3]]
## [[3]]$alpha
##  [1] 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548
##  [9] 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548
## [17] 47.6548 47.6548 47.6548 47.6548
## 
## [[3]]$beta
##  [1] 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428
##  [8] 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428
## [15] 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428
## 
## [[3]]$mu_a
## [1] 47.6548
## 
## [[3]]$mu_b
## [1] 5.239428
## 
## [[3]]$sigma
## [1] 15.41029  3.00000  3.00000

Pars <- c('mu_a', 'mu_b', 'sigma')

n.adapt <- 5000
n.update <-100000
n.iter <- 15000

jm = jags.model(data = M_data, "Meth_1.jags", n.chains = 3, inits = F_Inits, n.adapt = n.adapt)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 200
##    Unobserved stochastic nodes: 45
##    Total graph size: 867
## 
## Initializing model

update(jm, n.iter = n.update)
MCMCsamples = coda.samples(jm, n.iter = n.iter, variable.names = Pars, thin = 10)

## 
## Iterations = 105010:120000
## Thinning interval = 10 
## Number of chains = 3 
## Sample size per chain = 1500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##            Mean     SD Naive SE Time-series SE
## mu_a     37.646 1.7558 0.026174       0.166744
## mu_b     -4.747 3.0824 0.045950       0.281436
## sigma[1] 12.537 0.6338 0.009448       0.009591
## sigma[2]  0.815 0.6594 0.009829       0.029007
## sigma[3]  1.497 1.2381 0.018457       0.056728
## 
## 2. Quantiles for each variable:
## 
##               2.5%     25%     50%    75%  97.5%
## mu_a      34.37013 36.4124 37.6231 38.813 41.252
## mu_b     -11.00699 -6.8617 -4.6610 -2.629  1.123
## sigma[1]  11.35649 12.0906 12.5229 12.945 13.840
## sigma[2]   0.01857  0.2937  0.6685  1.176  2.436
## sigma[3]   0.05284  0.5576  1.2115  2.134  4.358

## Potential scale reduction factors:
## 
##          Point est. Upper C.I.
## mu_a           1.02       1.06
## mu_b           1.02       1.05
## sigma[1]       1.00       1.00
## sigma[2]       1.00       1.01
## sigma[3]       1.00       1.00
## 
## Multivariate psrf
## 
## 1.02

V. Parallel with different initial values

Repeat this in parallel. We use the worker PID to assign each worker one of the three sets of initial values. The results are similar to the sequential run, but let’s dig a little deeper.

cl <- parallel::makeCluster(3)

pid <- NA

for(i in 1:3){

  pidNum <- capture.output(cl[[i]])
  start <- regexpr("pid", pidNum)[[1]]
  end <- nchar(pidNum)
  pid[i] <- substr(pidNum, (start + 4), end)

}

parallel::clusterExport(cl, c("M_data", "n.adapt", "n.update", "n.iter", "Pars", "pid", "F_Inits"))

out <- clusterEvalQ(cl, {

  require(rjags)
  
  processNum <- which(pid==Sys.getpid())
  m_inits <- F_Inits[[processNum]]
  jm = jags.model(data = M_data, file = "Meth_1.jags", n.chains = 1, inits = m_inits, n.adapt = n.adapt)
  update(jm, n.iter = n.update)
  MCMCsamples = coda.samples(jm, n.iter = n.iter, variable.names = Pars, thin = 10)

  return(c(MCMCsamples, m_inits))

})

MCMCsamples <- mcmc.list(c(as.mcmc(out[[1]][1]), as.mcmc(out[[2]][1]), as.mcmc(out[[3]][1])))
(inits <- as.list(c(out[[1]][c(2:6)], out[[2]][c(2:6)], out[[3]][c(2:6)])))

## $alpha
##  [1] 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548
##  [9] 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548
## [17] 37.6548 37.6548 37.6548 37.6548
## 
## $beta
##  [1] -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572
##  [8] -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572
## [15] -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572
## 
## $mu_a
## [1] 37.6548
## 
## $mu_b
## [1] -4.760572
## 
## $sigma
## [1] 13.41029  1.00000  1.00000
## 
## $alpha
##  [1] 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548
##  [9] 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548
## [17] 42.6548 42.6548 42.6548 42.6548
## 
## $beta
##  [1] 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279
##  [8] 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279
## [15] 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279
## 
## $mu_a
## [1] 42.6548
## 
## $mu_b
## [1] 0.2394279
## 
## $sigma
## [1] 14.41029  2.00000  2.00000
## 
## $alpha
##  [1] 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548
##  [9] 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548
## [17] 47.6548 47.6548 47.6548 47.6548
## 
## $beta
##  [1] 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428
##  [8] 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428
## [15] 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428
## 
## $mu_a
## [1] 47.6548
## 
## $mu_b
## [1] 5.239428
## 
## $sigma
## [1] 15.41029  3.00000  3.00000

stopCluster(cl)

## 
## Iterations = 105010:120000
## Thinning interval = 10 
## Number of chains = 3 
## Sample size per chain = 1500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##             Mean     SD Naive SE Time-series SE
## mu_a     37.4562 2.0714 0.030879       0.205830
## mu_b     -4.5130 3.7122 0.055338       0.367036
## sigma[1] 12.5527 0.6529 0.009732       0.009559
## sigma[2]  0.8243 0.6413 0.009560       0.026329
## sigma[3]  1.4510 1.1574 0.017253       0.045274
## 
## 2. Quantiles for each variable:
## 
##               2.5%     25%     50%    75%  97.5%
## mu_a      33.60857 35.9598 37.5140 38.909 41.376
## mu_b     -11.71275 -7.1376 -4.5822 -1.935  2.614
## sigma[1]  11.36042 12.0948 12.5192 12.983 13.920
## sigma[2]   0.04813  0.3328  0.6741  1.170  2.374
## sigma[3]   0.04712  0.5425  1.1865  2.082  4.283

## Potential scale reduction factors:
## 
##          Point est. Upper C.I.
## mu_a           1.03       1.08
## mu_b           1.03       1.09
## sigma[1]       1.00       1.00
## sigma[2]       1.00       1.01
## sigma[3]       1.01       1.05
## 
## Multivariate psrf
## 
## 1.03

VII. Sequential with identical initial values

Let’s refit this model with identical initial values for each chain. We would never do this normally, but I want to illustrate a point. The chains are not identical, despite identical initial values. This is because .Random.seed will keep refreshing with new “random” values as JAGS works its way through the chains one by one.

set.seed(1) 
randomInits <- function( ){
  
  return(list(
    alpha = rnorm(SITES, 35, 1),
    beta = rnorm(SITES, 1, 1),
    mu_a = rnorm(1, 35, 2), 
    mu_b = rnorm(1, 35, 2), 
    sigma = runif(3, 1, 2)))
}

inits1 <- randomInits()
(F_Inits <- list(inits1, inits1, inits1))

## [[1]]
## [[1]]$alpha
##  [1] 34.37355 35.18364 34.16437 36.59528 35.32951 34.17953 35.48743
##  [8] 35.73832 35.57578 34.69461 36.51178 35.38984 34.37876 32.78530
## [15] 36.12493 34.95507 34.98381 35.94384 35.82122 35.59390
## 
## [[1]]$beta
##  [1]  1.9189774  1.7821363  1.0745650 -0.9893517  1.6198257  0.9438713
##  [7]  0.8442045 -0.4707524  0.5218499  1.4179416  2.3586796  0.8972123
## [13]  1.3876716  0.9461950 -0.3770596  0.5850054  0.6057100  0.9406866
## [19]  2.1000254  1.7631757
## 
## [[1]]$mu_a
## [1] 34.67095
## 
## [[1]]$mu_b
## [1] 34.49328
## 
## [[1]]$sigma
## [1] 1.757087 1.202692 1.711121
## 
## 
## [[2]]
## [[2]]$alpha
##  [1] 34.37355 35.18364 34.16437 36.59528 35.32951 34.17953 35.48743
##  [8] 35.73832 35.57578 34.69461 36.51178 35.38984 34.37876 32.78530
## [15] 36.12493 34.95507 34.98381 35.94384 35.82122 35.59390
## 
## [[2]]$beta
##  [1]  1.9189774  1.7821363  1.0745650 -0.9893517  1.6198257  0.9438713
##  [7]  0.8442045 -0.4707524  0.5218499  1.4179416  2.3586796  0.8972123
## [13]  1.3876716  0.9461950 -0.3770596  0.5850054  0.6057100  0.9406866
## [19]  2.1000254  1.7631757
## 
## [[2]]$mu_a
## [1] 34.67095
## 
## [[2]]$mu_b
## [1] 34.49328
## 
## [[2]]$sigma
## [1] 1.757087 1.202692 1.711121
## 
## 
## [[3]]
## [[3]]$alpha
##  [1] 34.37355 35.18364 34.16437 36.59528 35.32951 34.17953 35.48743
##  [8] 35.73832 35.57578 34.69461 36.51178 35.38984 34.37876 32.78530
## [15] 36.12493 34.95507 34.98381 35.94384 35.82122 35.59390
## 
## [[3]]$beta
##  [1]  1.9189774  1.7821363  1.0745650 -0.9893517  1.6198257  0.9438713
##  [7]  0.8442045 -0.4707524  0.5218499  1.4179416  2.3586796  0.8972123
## [13]  1.3876716  0.9461950 -0.3770596  0.5850054  0.6057100  0.9406866
## [19]  2.1000254  1.7631757
## 
## [[3]]$mu_a
## [1] 34.67095
## 
## [[3]]$mu_b
## [1] 34.49328
## 
## [[3]]$sigma
## [1] 1.757087 1.202692 1.711121

jm = jags.model(data = M_data, "Meth_1.jags", n.chains = 3, inits = F_Inits, n.adapt = n.adapt)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 200
##    Unobserved stochastic nodes: 45
##    Total graph size: 867
## 
## Initializing model

update(jm, n.iter = n.update)
MCMCsamples = coda.samples(jm, n.iter = n.iter, variable.names = Pars, thin = 10)

## 
## Iterations = 105010:120000
## Thinning interval = 10 
## Number of chains = 3 
## Sample size per chain = 1500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##             Mean     SD Naive SE Time-series SE
## mu_a     37.4203 1.7656 0.026320        0.15959
## mu_b     -4.4072 3.1421 0.046840        0.27421
## sigma[1] 12.5474 0.6371 0.009498        0.01031
## sigma[2]  0.8187 0.6232 0.009291        0.02418
## sigma[3]  1.5359 1.2003 0.017894        0.04630
## 
## 2. Quantiles for each variable:
## 
##               2.5%     25%     50%    75%  97.5%
## mu_a      33.90567 36.1916 37.4251 38.630 40.742
## mu_b     -10.34028 -6.5998 -4.4017 -2.281  1.761
## sigma[1]  11.34834 12.1032 12.5405 12.954 13.836
## sigma[2]   0.04406  0.3374  0.6834  1.165  2.334
## sigma[3]   0.08571  0.6393  1.2515  2.146  4.520

VI. Parallel with identical initial values

We now supply the same identical initial values to each of the R workers. This time however the resulting chains are identical! We will see why this is happening in the next section.

cl <- parallel::makeCluster(3)

pid <- NA

for(i in 1:3){

  pidNum <- capture.output(cl[[i]])
  start <- regexpr("pid", pidNum)[[1]]
  end <- nchar(pidNum)
  pid[i] <- substr(pidNum, (start + 4), end)

}

parallel::clusterExport(cl, c("M_data", "n.adapt", "n.update", "n.iter", "Pars", "pid", "randomInits", "SITES"))

out <- clusterEvalQ(cl, {

  require(rjags)
  
  processNum <- which(pid==Sys.getpid())
  set.seed(1)
  m_inits <- randomInits()
  jm = jags.model(data = M_data, file = "Meth_1.jags", n.chains = 1, inits = m_inits, n.adapt = n.adapt)
  update(jm, n.iter = n.update)
  MCMCsamples= coda.samples(jm, n.iter = n.iter, variable.names = Pars, thin = 10)

  return(c(MCMCsamples, m_inits))

})

MCMCsamples <- mcmc.list(c(as.mcmc(out[[1]][1]), as.mcmc(out[[2]][1]), as.mcmc(out[[3]][1])))
(inits <- as.list(c(out[[1]][c(2:6)], out[[2]][c(2:6)], out[[3]][c(2:6)])))

## $alpha
##  [1] 34.37355 35.18364 34.16437 36.59528 35.32951 34.17953 35.48743
##  [8] 35.73832 35.57578 34.69461 36.51178 35.38984 34.37876 32.78530
## [15] 36.12493 34.95507 34.98381 35.94384 35.82122 35.59390
## 
## $beta
##  [1]  1.9189774  1.7821363  1.0745650 -0.9893517  1.6198257  0.9438713
##  [7]  0.8442045 -0.4707524  0.5218499  1.4179416  2.3586796  0.8972123
## [13]  1.3876716  0.9461950 -0.3770596  0.5850054  0.6057100  0.9406866
## [19]  2.1000254  1.7631757
## 
## $mu_a
## [1] 34.67095
## 
## $mu_b
## [1] 34.49328
## 
## $sigma
## [1] 1.757087 1.202692 1.711121
## 
## $alpha
##  [1] 34.37355 35.18364 34.16437 36.59528 35.32951 34.17953 35.48743
##  [8] 35.73832 35.57578 34.69461 36.51178 35.38984 34.37876 32.78530
## [15] 36.12493 34.95507 34.98381 35.94384 35.82122 35.59390
## 
## $beta
##  [1]  1.9189774  1.7821363  1.0745650 -0.9893517  1.6198257  0.9438713
##  [7]  0.8442045 -0.4707524  0.5218499  1.4179416  2.3586796  0.8972123
## [13]  1.3876716  0.9461950 -0.3770596  0.5850054  0.6057100  0.9406866
## [19]  2.1000254  1.7631757
## 
## $mu_a
## [1] 34.67095
## 
## $mu_b
## [1] 34.49328
## 
## $sigma
## [1] 1.757087 1.202692 1.711121
## 
## $alpha
##  [1] 34.37355 35.18364 34.16437 36.59528 35.32951 34.17953 35.48743
##  [8] 35.73832 35.57578 34.69461 36.51178 35.38984 34.37876 32.78530
## [15] 36.12493 34.95507 34.98381 35.94384 35.82122 35.59390
## 
## $beta
##  [1]  1.9189774  1.7821363  1.0745650 -0.9893517  1.6198257  0.9438713
##  [7]  0.8442045 -0.4707524  0.5218499  1.4179416  2.3586796  0.8972123
## [13]  1.3876716  0.9461950 -0.3770596  0.5850054  0.6057100  0.9406866
## [19]  2.1000254  1.7631757
## 
## $mu_a
## [1] 34.67095
## 
## $mu_b
## [1] 34.49328
## 
## $sigma
## [1] 1.757087 1.202692 1.711121

stopCluster(cl)

## 
## Iterations = 105010:120000
## Thinning interval = 10 
## Number of chains = 3 
## Sample size per chain = 1500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##             Mean     SD Naive SE Time-series SE
## mu_a     37.6379 1.8668 0.027829       0.215473
## mu_b     -4.8200 3.5622 0.053102       0.441101
## sigma[1] 12.5690 0.6248 0.009315       0.009317
## sigma[2]  0.7183 0.5972 0.008903       0.025743
## sigma[3]  1.5330 1.1293 0.016835       0.042268
## 
## 2. Quantiles for each variable:
## 
##               2.5%     25%     50%    75%  97.5%
## mu_a      33.73693 36.4915 37.6967 39.007 41.104
## mu_b     -11.32854 -7.2022 -4.8806 -2.531  2.868
## sigma[1]  11.41017 12.1505 12.5393 12.983 13.858
## sigma[2]   0.01902  0.2192  0.5879  1.048  2.197
## sigma[3]   0.06657  0.6521  1.3063  2.175  4.114

VIII. Parallel with different initial values… is this enough?

Even if we make sure that the initial values are different for each chain running in parallel, we need to consider how JAGS handles randomization. JAGS sets its own seed (based on the time the model was implemented) and RNG (defaults to Mersenne-Twister) regardless of the RNG and seed set for the R worker it is running within. If each worker implements the model at the same time then each chain will have the same .Random.seed vector when performing MCMC. This is why the chains with identical initial values were also identical when run in parallel, even if the R worker seeds varied. If you have different initial values (or at least one different initial initial value for a free parameter in your model) then the chains will not be identical due to the MCMC algorithm itself. You can specify a seed and RNG for each chain directly with its initial values. For example, in addition to different initial values, we specify different starting seeds and RNGs for each chain.

inits1 <- list(
  alpha = coef(M1)$sites.lmer[,1], 
  beta = coef(M1)$sites.lmer[,2],
  mu_a = summary(M1)[[10]][1, 1], 
  mu_b = summary(M1)[[10]][2, 1], 
  sigma = sigma + 1,
  .RNG.name = "base::Mersenne-Twister",
  .RNG.seed = 1)
    
inits2 <- list(
  alpha = coef(M1)$sites.lmer[,1] + 5, 
  beta = coef(M1)$sites.lmer[,2] + 5,
  mu_a = summary(M1)[[10]][1, 1] + 5, 
  mu_b = summary(M1)[[10]][2, 1] + 5, 
  sigma = sigma + 2,
  .RNG.name = "base::Marsaglia-Multicarry",
  .RNG.seed = 22)

inits3 <- list(
  alpha = coef(M1)$sites.lmer[,1] + 10, 
  beta = coef(M1)$sites.lmer[,2] + 10,
  mu_a = summary(M1)[[10]][1, 1] + 10, 
  mu_b = summary(M1)[[10]][2, 1] + 10, 
  sigma = sigma + 3,
  .RNG.name = "base::Wichmann-Hill",
  .RNG.seed = 373)

F_Inits <- list(inits1, inits2, inits3)

cl <- parallel::makeCluster(3)

pid <- NA

for(i in 1:3){

  pidNum <- capture.output(cl[[i]])
  start <- regexpr("pid", pidNum)[[1]]
  end <- nchar(pidNum)
  pid[i] <- substr(pidNum, (start + 4), end)

}

parallel::clusterExport(cl, c("M_data", "n.adapt", "n.update", "n.iter", "Pars", "pid", "F_Inits"))

out <- clusterEvalQ(cl, {

  require(rjags)
  
  processNum <- which(pid==Sys.getpid())
  m_inits <- F_Inits[[processNum]]
  jm = jags.model(data = M_data, file = "Meth_1.jags", n.chains = 1, inits = m_inits, n.adapt = n.adapt)
  update(jm, n.iter = n.update)
  MCMCsamples = coda.samples(jm, n.iter = n.iter, variable.names = Pars, thin = 10)

  return(c(MCMCsamples, m_inits))

})

MCMCsamples <- mcmc.list(c(as.mcmc(out[[1]][1]), as.mcmc(out[[2]][1]), as.mcmc(out[[3]][1])))
(inits <- as.list(c(out[[1]][c(2:6)], out[[2]][c(2:6)], out[[3]][c(2:6)])))

## $alpha
##  [1] 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548
##  [9] 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548 37.6548
## [17] 37.6548 37.6548 37.6548 37.6548
## 
## $beta
##  [1] -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572
##  [8] -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572
## [15] -4.760572 -4.760572 -4.760572 -4.760572 -4.760572 -4.760572
## 
## $mu_a
## [1] 37.6548
## 
## $mu_b
## [1] -4.760572
## 
## $sigma
## [1] 13.41029  1.00000  1.00000
## 
## $alpha
##  [1] 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548
##  [9] 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548 42.6548
## [17] 42.6548 42.6548 42.6548 42.6548
## 
## $beta
##  [1] 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279
##  [8] 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279
## [15] 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279 0.2394279
## 
## $mu_a
## [1] 42.6548
## 
## $mu_b
## [1] 0.2394279
## 
## $sigma
## [1] 14.41029  2.00000  2.00000
## 
## $alpha
##  [1] 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548
##  [9] 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548 47.6548
## [17] 47.6548 47.6548 47.6548 47.6548
## 
## $beta
##  [1] 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428
##  [8] 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428
## [15] 5.239428 5.239428 5.239428 5.239428 5.239428 5.239428
## 
## $mu_a
## [1] 47.6548
## 
## $mu_b
## [1] 5.239428
## 
## $sigma
## [1] 15.41029  3.00000  3.00000

stopCluster(cl)

## 
## Iterations = 105010:120000
## Thinning interval = 10 
## Number of chains = 3 
## Sample size per chain = 1500 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##             Mean     SD Naive SE Time-series SE
## mu_a     37.4929 1.7477 0.026053       0.145026
## mu_b     -4.5576 3.1353 0.046738       0.285539
## sigma[1] 12.5331 0.6426 0.009579       0.009581
## sigma[2]  0.8366 0.6593 0.009828       0.027381
## sigma[3]  1.5739 1.2601 0.018784       0.058596
## 
## 2. Quantiles for each variable:
## 
##               2.5%     25%     50%    75%  97.5%
## mu_a      34.18909 36.2443 37.4871 38.768 40.726
## mu_b     -10.37801 -6.8363 -4.5603 -2.357  1.459
## sigma[1]  11.33913 12.0772 12.5109 12.965 13.879
## sigma[2]   0.02438  0.3218  0.6996  1.197  2.474
## sigma[3]   0.04055  0.5762  1.2948  2.253  4.720

## Potential scale reduction factors:
## 
##          Point est. Upper C.I.
## mu_a           1.01       1.02
## mu_b           1.01       1.02
## sigma[1]       1.00       1.00
## sigma[2]       1.00       1.01
## sigma[3]       1.01       1.02
## 
## Multivariate psrf
## 
## 1.01